
# Make Figure 3, and calculate the uncertainty ranges reported in column 6 of Table 2. 

library(class)  #user will need to install class package if not already installed

wd <- ""  #insert into quotation marks the working directory to location where replication file was unzipped. E.g /Documents/BDLMS_Replication/ 	Make sure to include the final backslash.
"%&%"<-function(x,y)paste(x,y,sep="")  #define a function for easy string pasting


numstud=7
strip=vector("list",numstud)  #set number of different studies as second value here. this is going to be a list that we fill with bootstrapped future projections for each study. then we loop over the list to make the Figure. 

#####################################
# MNS and SHF
#####################################

setwd(wd%&%"data/MNS_SHF/")  #navigate to MNS and SHF data directory
mns=read.csv("boot_mns.csv")  #bootstrapped regression coefficients, column 2 in Table 3 of MNS
shf=read.csv("boot_shf.csv")  # same specification, but restricted to non-irrigated non-urban counties as in SHF
data=read.csv("mns_sample.csv")
load("MNS_SHF_climchg_temp.Rdata")  #temperature projections for MNS and SHF
tm=tm*9/5  # MNS and SHF temperature variables are in farenheight, but climate change projections are in C
load("MNS_SHF_climchg_prec.Rdata")  #precip projections for MNS and SHF
baseprec=data.frame(data$janprec_raw,data$aprprec_raw,data$julprec_raw,data$octprec_raw)

p=1 #set period at mid-century.

# MNS first
inlist=7  #this is where to put it in the list that we are filling
smpl=(is.na(data$cropsales)==F)
m1=as.matrix(mns[,2:17])  #boostrapped estimates for the 16 independent variables - in $/acre
m2=c()  #the values at which to evaluate the boootstrapped estimates
for (i in 1:4) {
	m2=rbind(m2,tm[,i,p,1],tm[,i,p,1]^2)
	}
for (i in 1:4) {
	prchg=pr[,i,p,1]*weighted.mean(baseprec[smpl,i],data$farmland[smpl])  #projections are in %chg, converting to levels to evaluate coeffs
	m2=rbind(m2,prchg,prchg^2)
	}

medboot=apply(m1,2,median)  #median bootstrapped coefficients, about the same as reg coefficients
imp=medboot%*%m2/weighted.mean(data$farmval1982[smpl],data$cropsales[smpl])*100  #impact projection for each model, converted to % change, and truncated at -100% losses
imp=replace(imp,imp<(-100),-100)
mods=dimnames(tm)[[1]]  # names of the different climate models
# now we are going to fill the list element with some subelements that we want:  the projections for our 3 scenarios, the name of the study, its year, and the outcome it looks at.  we do the same thing for all the other studies below.
strip[[inlist]][[1]]=imp[grep("a1b",mods)]  
strip[[inlist]][[2]]=imp[grep("a2",mods)]
strip[[inlist]][[3]]=imp[grep("b1",mods)]
strip[[inlist]][[4]]="MNS"
strip[[inlist]][[5]]=1994
strip[[inlist]][[6]]="land val."


# SHF
inlist=6
smpl=data$dryrural==1  #the sample of non-irrigated rural counties we are using
m1=as.matrix(shf[,2:17])  #boostrapped estimates for the 16 independent variables - in $/acre
m2=c()  #the values at which to evaluate the boootstrapped estimates
for (i in 1:4) {
	m2=rbind(m2,tm[,i,p,2],tm[,i,p,2]^2)
	}
for (i in 1:4) {
	prchg=pr[,i,p,2]*weighted.mean(baseprec[smpl,i],data$farmland[smpl])  #projections are in %chg, converting to levels to evaluate coeffs
	m2=rbind(m2,prchg,prchg^2)
	}

medboot=apply(m1,2,median)  #median bootstrapped coefficients, about the same as reg coefficients
imp=medboot%*%m2/weighted.mean(data$farmval1982[smpl],data$totalcropland1982[smpl])*100 #impact projection for each model, converted to % change, and truncated at -100% losses
imp=replace(imp,imp<(-100),-100)
mods=dimnames(tm)[[1]]  #modelnames
strip[[inlist]][[1]]=imp[grep("a1b",mods)]
strip[[inlist]][[2]]=imp[grep("a2",mods)]
strip[[inlist]][[3]]=imp[grep("b1",mods)]
strip[[inlist]][[4]]="SHF"
strip[[inlist]][[5]]=2005
strip[[inlist]][[6]]="land val."


##############################################################
# 	DJO
##############################################################

inlist=1

setwd(wd%&%"data/DJO/")  #navigate to DJO data directory

load("DJO_projections.Rdata")  #the impact projections we calculated
boot=read.csv("boot_djo.csv")  #boostrapped reg coefficients from DJO specification in Table 2, Col 3
boot=boot[,2] #just keep temp coefficients

#collect scenario and model names
nms=strsplit(dimnames(impact)[[2]],"\\_")
scen=mod=per=c()
for (i in 1:length(nms)) {
	scen=c(scen,nms[[i]][1]) #scenario
	mod=c(mod,nms[[i]][2]) #climate model
	per=c(per,nms[[i]][3]) #mid- or end-of-century
}

# location of median bootstrap temperature coefficient in impact array
medloc=as.numeric(knn1(boot,median(boot),1:1000))
med=impact[medloc,]

strip[[inlist]][[1]]=med[per=="mid"&scen=="a1b"]
strip[[inlist]][[2]]=med[per=="mid"&scen=="a2"]
strip[[inlist]][[3]]=med[per=="mid"&scen=="b1"]
strip[[inlist]][[4]]="DJO"
strip[[inlist]][[5]]=2012
strip[[inlist]][[6]]="GDP/cap"


##############################################################
# 	DG 2007
##############################################################

setwd(wd%&%"data/DG/")  #navigate to DG directory

inlist=5

load("DG_climchg.Rdata")
boot=read.csv("boot_dg.csv")
data=read.csv("dg_data.csv")

m1=as.matrix(boot[,2:9])
tot=aggregate(data,by=list(data$fips),mean,na.rm=T)
dryacre=sum(tot$fland[tot$dry==1])
irracre=sum(tot$fland[tot$dry==0])
totprof=sum(tot$fland*tot$y,na.rm=T)

#coefficients are in $/acre, so dividing by average dryland or irrigated profits/acre to get impacts in %.
m1[,1:4]=m1[,1:4]*dryacre/totprof*100
m1[,5:8]=m1[,5:8]*irracre/totprof*100

mods=dimnames(climchg)[[1]]  #modelnames
medboot=apply(m1,2,median)  #median bootstrapped coefficients, about the same as reg coefficients
m2=t(climchg[,,1])  #mid century climate change projections
med=medboot%*%m2  #climate uncertainty, evaluated at median bootstrapped coefficients

strip[[inlist]][[1]]=med[grep("a1b",mods)]
strip[[inlist]][[2]]=med[grep("a2",mods)]
strip[[inlist]][[3]]=med[grep("b1",mods)]
strip[[inlist]][[4]]="DG"
strip[[inlist]][[5]]=2007
strip[[inlist]][[6]]="ag profits"


##############################################################
# 	FHRS 2010
##############################################################

setwd(wd%&%"data/FHRS/")  #navigate to FHRS directory

inlist=2

load("FHRS_climchg.Rdata")
boot=read.csv("boot_fhrs.csv")
data=read.csv("fhrs_data.csv")

tot=aggregate(data,by=list(data$fips),mean,na.rm=T)
totprod=sum(tot$corn_planted*tot$yield,na.rm=T)

# multiply area-weighted yield effect (what the regression spits out) by corn area to get effect of change on overall production, then divide by total production
m1=as.matrix(boot[,2:9])
m1[,1:4]=m1[,1:4]*sum(tot$corn_planted[tot$dry==1],na.rm=T)/totprod*100
m1[,5:8]=m1[,5:8]*sum(tot$corn_planted[tot$dry==0],na.rm=T)/totprod*100

mods=dimnames(climchg)[[1]]  #modelnames
medboot=apply(m1,2,median)  #median bootstrapped coefficients
m2=t(climchg[,,1])  #mid century climate change projections
med=medboot%*%m2  #climate uncertainty, evaluated at median bootstrapped coefficients

strip[[inlist]][[1]]=med[grep("a1b",mods)]
strip[[inlist]][[2]]=med[grep("a2",mods)]
strip[[inlist]][[3]]=med[grep("b1",mods)]
strip[[inlist]][[4]]="FHRS"
strip[[inlist]][[5]]=2010
strip[[inlist]][[6]]="corn yield"


##############################################################
# 	BMSDL 2009
##############################################################

inlist=4
setwd(wd%&%"data/BMSDL/")  #navigate to BMSDL directory

boot=read.csv("boot_BMSDL.csv")  #boostrapped reg coefficients from pnas specification 1 in Table 1 of PNAS paper
temp=read.csv("BMSDL_climchg.csv")  #temperature projections

# convert coefficients to %change
boot=boot/0.11*100  #0.11 is the baseline incidence of conflict in the sample
tchg=apply(temp[,4:dim(temp)[2]],2,mean)

#collect scenario and model names
nms=strsplit(names(tchg),"\\_")
scen=mod=per=c()
for (i in 1:length(nms)) {
	scen=c(scen,nms[[i]][1]) #scenario
	mod=c(mod,nms[[i]][2]) #climate model
	per=c(per,nms[[i]][3]) #mid- or end-of-century
}

# median bootstrap coefficients
medboot = apply(boot[,2:3],2,median)
m2=as.matrix(t(data.frame(tchg,tchg)),ncol=length(tchg))  ##two temp coefficients that we're evaluating together
med=medboot%*%m2  #climate uncertainty

strip[[inlist]][[1]]=med[per=="mid"&scen=="a1b"]
strip[[inlist]][[2]]=med[per=="mid"&scen=="a2"]
strip[[inlist]][[3]]=med[per=="mid"&scen=="b1"]
strip[[inlist]][[4]]="BMSDL"
strip[[inlist]][[5]]=2009
strip[[inlist]][[6]]="civil war"


##############################################################
# 	SL 2010
##############################################################

inlist=3
setwd(wd%&%"data/SL/")  #navigate to SL directory

clim=read.csv('SL_climchg.csv')
boot=read.csv('boot_sl.csv')
data=read.csv('africa_yield_clim.csv')

prec_base=aggregate(data$maize_prec,by=list(data$fidcode),mean,na.rm=T)  #country-average precip
prec_base=mean(prec_base[,2])

tchg=apply(clim[,2:103],2,mean)  #regressions are unweighted so evaluating at a simple average of temp and precip change across countries
pchg=apply(clim[,104:dim(clim)[2]],2,mean)*prec_base  #precip changes are in %, so evaluating at continental mean annual precip to get changes in levels. 

#collect scenario and model names
nms=strsplit(names(tchg),"\\_")
scen=mod=per=c()
for (i in 1:length(nms)) {
	scen=c(scen,nms[[i]][2]) #scenario
	mod=c(mod,nms[[i]][3]) #climate model
	per=c(per,nms[[i]][4]) #mid- or end-of-century
}

m1=data.matrix(boot[,2:3])  #bootstrapped coefficients
medboot = apply(m1,2,median)
m2=data.matrix(rbind(tchg,pchg))  #clim chg projections
med=medboot%*%m2*100

strip[[inlist]][[1]]=med[per=="mid"&scen=="a1b"]
strip[[inlist]][[2]]=med[per=="mid"&scen=="a2"]
strip[[inlist]][[3]]=med[per=="mid"&scen=="b1"]
strip[[inlist]][[4]]="SL"
strip[[inlist]][[5]]=2009
strip[[inlist]][[6]]="corn yield"


##############################################################
# 	Now make Figure 3
##############################################################

setwd(wd%&%"output/")  #navigate to SL directory
sc=c("A1B","A2","B1")
#pdf(file="Figure3.pdf",width=6,height=8)
par(mar=c(4,7,2,1))
plot(1,xlim=c(-120,120),ylim=c(0.7,numstud+0.3),type="n",yaxt="n",xlab="% impact",ylab="",xaxt="n")
sq=seq(-100,125,25)
axis(1,at=sq,sq,cex.axis=0.8)
abline(v=0,lwd=1.5,lty=2)
abline(v=sq,lty=3,col="grey")
abline(h=1:(numstud-1)+0.5)
ord=c(3,1,2)
for (i in 1:numstud) {
	for (j in 1:3) {
		d=0.25
		loc=i+seq(-d,d,d)[ord[j]]
		pts=strip[[i]][[j]]
		points(pts,rep(loc,length(pts)),cex=1.3,col="grey",pch=124)
		if (j==1) {
			points(pts[18],loc,cex=1.5,col="black",pch=124)  #plot a1b hadley as darker line
			}
		text(-115,loc,sc[j],cex=0.7)
		}
	mtext(strip[[i]][[4]],2,las=1,line=1,cex=0.8,at=i+0.1)
	mtext(strip[[i]][[6]],2,las=1,line=1,cex=0.75,at=i-0.05)
}
#dev.off()
dev.copy2eps(file="Figure3.eps",width=6,height=8)

##############################################################
# 	Calculate ranges for table 2
##############################################################

for (j in 1:numstud) {
	clim=round(unlist(strip[[j]][1]),0)
	print(paste("Range in ",strip[[j]][4]," = ",min(clim),",",max(clim),sep=""))
}
